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ABSTRACT 

We address complexity issues for linear differential equa- 
tions in characteristic p > 0: resolution and computation of 
the p-curvature. For these tasks, our main focus is on al- 
gorithms whose complexity behaves well with respect to p. 
We prove bounds linear in p on the degree of polynomial 
solutions and propose algorithms for testing the existence of 
polynomial solutions in sublinear time 0(p 1,/2 ), and for de- 
termining a whole basis of the solution space in quasi-linear 
time 0(p); the O notation indicates that we hide logarithmic 
factors. We show that for equations of arbitrary order, the 
p-curvature can be computed in subquadratic time 0(p 179 ), 
and that this can be improved to 0(log(p)) for first order 
equations and to 0(p) for classes of second order equations. 

Categories and Subject Descriptors: 
1.1.2 [Computing Methodologies]: Symbolic and Alge- 
braic Manipulation - Algebraic Algorithms 

General Terms: Algorithms, Theory 

Keywords: Algorithms, complexity, differential equations, 
polynomial solutions, p-curvature. 

1. INTRODUCTION 

We study several algorithmic questions related to linear dif- 
ferential equations in characteristic p, where p is a prime 
number: resolution of such equations and computation of 
their p-curvature. Our emphasis is on the complexity view- 
point. 

Let thus ¥ p be the finite field with p elements, and let 
¥ p (x){d) be the algebra of differential operators with coeffi- 
cients in ¥ p (x), with the commutation relation dx = xd+ 1. 
One of the important objects associated to a differential op- 
erator L of order r in ¥ p (x){d) is its p-curvature, hereafter 
denoted A p . By definition, this is the (r x r) matrix with 
coefficients in ¥ p (x), whose (i,j)-entry is the coefficient of 
d % in the remainder of the Euclidean (right) division of d p+3 
by L, for < i,j < r. 
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The concept of p-curvature originates in Grothendieck's 
work in the late 1960s, in connection to one of his famous 
(still unsolved) conjectures. In its simplest form, this con- 
jecture is an arithmetic criterion of algebraicity, which states 
that a linear differential equation with coefficients in Q(x) 
has a basis of algebraic solutions over Q(x) if and only if 
its reductions modulo p have zero p-curvature, for almost 
all primes p. The search of a proof of this criterion moti- 
vated the development of a theory of differential equations 
in characteristic p by Katz [20] , Dwork [14] , Honda [19] , etc. 

There are two basic differences between differential equa- 
tions in characteristic zero and p: one concerns the dimen- 
sion of the solution space, the other, the form of the so- 
lutions. While in characteristic zero, a linear differential 
equation of order r admits exactly r linearly independent 
solutions, this is no longer true in positive characteristic: 
for L £ ¥ p (x){d), the dimension of the solution space of 
the equation Ly = over the field of constants ¥ p (x p ) is 
generally less than the order r. Moreover, by a theorem of 
Cartier and Katz (see Lemma 2 below), the dimension is 
exactly r if and only if the p-curvature matrix A p is zero. 
Thus, roughly speaking, the p-curvature measures to what 
extent the solution space of a differential equation modulo p 
has dimension close to its order. 

On the other hand, the form of the solutions is simpler 
in characteristic p than in characteristic zero. Precisely, the 
existence of polynomial solutions is equivalent to the exis- 
tence of solutions which are either algebraic over ¥ p (x), or 
power series in F p [[a;]], or rational functions in ¥ p (x) [19]. 
Therefore, in what follows, by solving Ly = we simply 
understand finding its polynomial solutions. 

In computer algebra, the p-curvature was publicised by 
van der Put [25, 26], who used it as a central tool in designing 
algorithms for factoring differential operators in ¥ p (x)(d). 
Recently, his algorithms were analyzed from the complexity 
perspective and implemented by Cluzeau [11], who extended 
them to the case of systems. Cluzeau also took in [12] a first 
step towards a systematic modular approach to the algorith- 
mic treatment of differential equations. 

Improving the complexity of the p-curvature computation 
is an interesting problem in its own right. Our main moti- 
vation for studying this question comes, however, from con- 
crete applications. First, in a combinatorial context, the use 
of the p-curvature served in the automatic classification of 
restricted lattice walks [8] and notably provided crucial help 
in the treatment of the notoriously difficult case of Gessel's 
walks [7]. Also, intensive p-curvature computations were 



needed in [4], where the question is to decide whether var- 
ious differential operators arising in statistical physics have 
nilpotent, or zero, p-curvature. 

In the latter questions, the prime p was "large" , typically 
of the order of 10 4 . This remark motivates our choice of 
considering p as the most important parameter: our primary 
objective is to obtain complexity estimates featuring a low 
exponent in p. 

Previous work. The non-commutativity of ¥ p (x){d) pre- 
vents one from straightforwardly using binary powering tech- 
niques for the computation of A p via that of d p mod L. 
Thus, the complexity of all currently known algorithms for 
computing the p-curvature is quadratic in p. 

Katz [21] gave the first algorithm, based on the following 
matrix recurrence: define 

Ai = A, A fc+1 = A' fc + AA fe , (1) 

where A € ^ r (¥ p (x)) is the companion matrix associated to 
L; then, A p is the p-curvature matrix (hence our notation). 

It was observed in [27, §13.2.2] that it is slightly more effi- 
cient to replace (1) by the recurrence Vfe+i = v' k +A.v k which 
computes the first column v/t of A&, by taking for vo the first 
column of I r . Then v p , . . . , v p + r _i are the columns of A p . 
This alternative requires only matrix-vector products, and 
thus saves a factor of r, but still remains quadratic in p. 
Cluzeau proposed in [11, Prop. 3.2] a fraction-free version 
of (1) having essentially the same complexity, but incor- 
rectly stated that the method in [27] works in linear time 
in p. 

Concerning polynomial and rational solutions of differen- 
tial equations modulo p, very few algorithms can be found in 
the literature. Cluzeau proposes in [11, §2] an algorithm of 
cubic complexity in p and, in the special case when A p = 0, 
a different algorithm of quadratic complexity in p, based on 
a formula due to Katz which is the nub of Lemma 2 below. 

Our contribution. We prove in Section 3 a linear bound 
in p on the degree for a basis of the solution space of poly- 
nomial solutions of an equation Ly — 0. Then, we adapt 
the algorithm in [1] and its improvements [6] to the case of 
positive characteristic; we show how to test the existence 
of polynomial solutions in time nearly proportional to p 1 ^ 2 , 
and how to determine a full basis of the solution space in 
time quasi-linear in p. 

Regarding the p-curvature, we first focus on two particular 
cases: first order operators, where the cost is polynomial 
in log(p) (Section 4), and second order ones, for which we 
obtain a cost quasi-linear in p in some cases (Section 5). 

In general, a useful way to see (1) is to note that the p- 
curvature is obtained by applying the operator (d + A) p_1 
to A. In Section 6 we exploit this observation. As a side 
result, we give a baby steps / giant steps algorithm for com- 
puting the image Lu of an operator L applied to a polyno- 
mial u; this is inspired by Brent-Kung's algorithm for power 
series composition [9], 

Complexity measures. Time complexities are measured 
in terms of arithmetic operations in F p . 

We let M : N — » N be such that polynomials of degree 
at most n in F p [a;] can be multiplied in time M(n). Fur- 
thermore, we assume that M(n) satisfies the usual assump- 
tions of [18, §8.3]; using Fast Fourier Transform, M(n) can 
be taken in 0(nlogn log log n) [23, 10]. We suppose that 
2 < cu < 3 is a constant such that two matrices in ^d(F p ) 



can be multiplied in time 0(n w ). The current tightest upper 
bound is u < 2.376 [13]. 

The precise complexity estimates of our algorithms are 
sometimes quite complex; to highly their main features, we 
rather give simplified estimates. Thus, we use the notation 
/ £ 6(g) for f,g : N -> N if / is in 0(glog(g) m ) for some 
m > 1. For instance, M(n) is in O(n). 

2. PRELIMINARIES 

Basic properties of the p-curvature. We first give de- 
gree bounds on the p-curvature of an operator. Consider 

L = i (x) + £ 1 {x)d + ■■■+£ r (x)d r , (2) 

with all £i in ¥ p [x] of degrees at most d and £ r ^ 0. As 
in (1), we define Ai = A and Afc+i = A' k + AAfc for k > 1. 

Lemma 1. For k > 0, let B k = i*A k . Then B k is in 
-# r (¥ v [x]), with entries of degree at most dk. 

Proof. Explicitly, we have 




From this, we see that the sequence Bfe satisfies the equation 

B fc+1 = e r B' h + (Bi - ke%)B k , 

where I r is the r x r identity matrix. The claim follows. □ 
In particular, the p-curvature A p has the form B p /££, with 
B p a polynomial matrix of degree at most dp. 

A second useful result is the following lemma, attributed 
to Katz. It relates the solution space of Ly — to the p- 
curvature and generalizes a theorem of Cartier. A proof can 
be found in [11, Th. 3.8]. 

Lemma 2. The dimension over¥ p (x p ) of the vector space 
of rational solutions of L is equal to the dimension over 
¥ p (x) of the kernel of A p . In particular, L has a basis of 
polynomial solutions if and only if its p-curvature is zero. 

Operator algebras. In what follows, we mainly consider 
operators with coefficients in F p (a;), but also sometimes more 
generally in the (n x n) matrix algebra .^ n (F p (x)); as has 
been done up to now, we will write matrices in bold face. If 
L is in ^ n (¥ p {x))(d) of the form 

L = £ (x) + £j (x)d +■■■ + £ r (x)d r , 

with coefficient matrices li in ^( n (¥ p [x]) of maximal degree 
d, we say that L has bidegree (d,r). 

Regularization. For most of our algorithms, we must as- 
sume that the origin x — does not cancel the leading term 
£ r £ ¥ p [x] of the operator L. 

If we can find x £ F p such that £ r {x') ^ 0, we can ensure 
this property by translating the origin to x' . To ensure 
that we can find x' , we must make the following hypothesis, 
written H: £ r does not vanish identically on F p . 

Lemma 3. Given L of bidegree (d,r), testing whether H 
holds can be done in time 0(M(d)) C O(d). If so, one can 
find x' such that £ r (x') 7^ and translate the coordinates' 
origin to x' in time 0(rM(d) log(d)) C 0(rd). 



Proof. Testing H amounts to verify whether x p — x divides 
l r . If deg(£ r ) < p, H obviously holds. Else, we have p < d; 
then, it is enough to reduce £ r modulo x p — x, which takes 
time 0(M{d)). 

If H holds, we know that we can find x' £ {0, . . . , deg(£ r )} 
such that £ r (x') 7^ 0; so it is enough to evaluate i r at this 
set of points, which by [18, §10.1] takes time 0(M(d) log(d)). 
Once x' is known, we shift all coefficients of L by x' . Us- 
ing fast algorithms for polynomial shift [17], the time is 
0(M(d) log(d)) per coefficient; the conclusion follows. □ 

As a consequence, in all the following algorithms, we will 
assume that H holds. If not, one could actually work in a 
low-degree extension of F p to find x'\ we do not consider this 
generalization here. 

3. POLYNOMIAL SOLUTIONS 

We start with the study of the polynomial solutions of a 
linear differential equation; aside from its own interest, this 
question will arise in our algorithm for order two operators 
in Section 5. 

Theorem 1. Let L be as in (2), with r < d and r < 
p, and such that H holds. Then, one can test whether the 
equation Lu = has non-zero solutions in F p (a;) in time 

,U1 1/2 1/2 . jul + l 

0(d r ' p ' + a r ). 

If so, one can determine a basis of the solution set consisting 
of polynomials of degree at most dp — 1 in extra time 

0(d u+1 rp + d 2 r u ' +3 p). 

The main point here is that for fixed d and r, testing the 
existence of solutions takes time Ofj} 1 ' 2 ), whereas finding a 
basis of the solution space takes time 0(p). 

In all this section, L is fixed, and the assumptions of The- 
orem 1 are satisfied. The assumptions on the relative order 
of magnitude of p, d, r help us obtain simple cost estimates 
and rule out some possible overlaps in indices modulo p. 
The assumption r < d is here mostly for convenience; the 
assumption r < p is necessary. 

3.1 Degree bounds 

Let T be the F p (i p )-vector space of rational solutions of 
the equation Lu — 0. The following proposition proves a 
bound linear in p on the degree of a basis of T. To our 
knowledge, such linear bounds were previously available only 
in two particular cases: (a) when the equation has a basis of 
polynomial solutions and under the additional hypotheses 
< deg(£ ) - r < p - 1 and p > r [19, Th. 7]; (b) when 
r = 2 and the equation has exactly one nonzero polynomial 
solution [14, Lemma 10.1]. These bounds are respectively 
(p— r)d+(X) for (a) and i(p— l)(d— 1) for (b). In the general 
case, the analysis in [11, 12] suggests a bound quadratic in 
p of type p(p + d). Our result refines this approach. 

Proposition 1. If Lu = has at least one nonzero solu- 
tion in¥ p (x), then J- admits a basis consisting of polynomial 
solutions of degree at most pd — 1 each. 

PROOF. The map ifiL ■ ¥ p (x) — + ¥ p (x) defined by y 1— » L(y) 
is ¥ p (x p )- linear. Let M £ ^£ p (¥ p (x p )) be the matrix of this 
map with respect to the basis (l,x, . . . ,x p ~ 1 ). Write M — 
(mij) <i,j< p -i for some rm,j in F p [a; p ]. Then, u £ ¥ p [x] is 



in T if and only if M x [«o • ■ ■ it p -i] f = 0, with Ui in F p [s p ] 
such that u — uo + u\x + ••• + u p -ix p ~ . 

Since L(x x ) = X/ < -1 m i,i x ° is a sum of p polynomials 
of pairwise distinct degrees deg(mij) + j, we deduce that 
for all i, j, deg(m 4j ) + j < deg(L(x 1 )). 

Since Lu — has a non-zero solution in F p (a;), it has 
also a non-zero solution in F p [s], by clearing denomina- 
tors. Let thus v be in ¥ p [x] \ {0} such that Lv = 0, or 
equivalently £ov = — X^i<j< r ^j'V ■ Since all terms in the 
right-hand side have degree at most d + deg(v) — 1, we de- 
duce that deg(£o) < d — 1. This implies that L(x z ) — 
4s'+Ei<j< r * " ' (i—j+l)£i xl ~ :> nas degree at most d+i— 1. 

To summarize, for all < i, j < p — 1, we obtain the 
inequality deg(mi J ) < (d — 1) + (i — j). This implies that 
for any permutation a of {0, . . . , p — 1}, 

de s(nr=o m i,T(i)) = Efco de g( m i,<r(o) < p( d - !)> 

since the sum of the terms i — a(i) is zero. This implies that 
all minors of M have degree at most p(d— 1), since any term 
appearing in the expansion of such minors can be completed 
to form one of the form Ilo<i<p-i m ;,<r(i)- 

The nullspace of M admits a basis [vi , . . . , v/t] , all of 
whose entries are minors of M. By what was said above, 
they all have degree at most p(d — 1). A basis of T is 
easily deduced: to v; = [v ^0 • ■ • Mi,p-i] corresponds the 
polynomial Vi = u^o + • • • + Vi, P -ix p ~ 1 . We deduce that 
deg(ui) < p — 1 + p(d — 1) = pd — 1, as claimed. □ 

3.2 Solutions of bounded degree 

Let Q C F p [a:] be the F p -vector space of polynomial solutions 
of Lu = of degree at most pd — 1. We are interested in 
computing either the dimension of Q, or an F p -basis of it. In 
view of the former proposition, this will be sufficient to prove 
Theorem 1. Proposition 2 gives cost estimates for these 
tasks, adapting the algorithm in [1] and its improvements [6]. 

Proposition 2. Under the assumptions of Theorem 1, 
one can compute dimf p (Q) in time 

/\/ 1UJ 1/2 1/2 1 JU+1 UJ-1\ 

(J(d r ' p ' + a r ). 

One can deduce a basis of Q in extra time 0(d" + rp). 

For r and d fixed, the main feature of this result is that 
the cost of computing the dimension of Q is the sublinear 
(5(p 1//2 ), whereas the cost of computing a basis of it is 0(p). 

Proof. Let uo,...,u p d-i be unknowns and let u be the 
polynomial u = ~^2 n<pd u n x n ; for n < or n > pd, we let 
u n = 0. There exist Co, • • • , Cd+ r in F p [n], of degree at most 
r, such that for n > 0, the coefficient of degree n of 

£ (x)u+---+£ r (x)u {r) (3) 

is C„ = co(n)u n -d + ■ ■ ■ + Cd+ r (n)u n + r ; note for further use 
that 

C„- r = c (n — r)u n _ r _d H h c d+r {n - r)u n . (4) 

The polynomial u is in Q if and only if C„ = for < n < 
(p + l)d — 1. Shifting indices, we obtain the system of linear 
equations C n - r = 0, with r < n < (p + l)d + r — 1, in the 
unknowns uo, ■ ■ ■ , u pd -\ . 

The matrix of this system is band-diagonal, with a band of 
width d+r+1. In characteristic zero or large enough, one can 



eliminate each unknown u n , with n > r, using C„_ r . Here, 
some equations C„_ r become deficient, in the sense that the 
coefficient of u n vanishes; this induces a few complications. 

Outline of the computation. Since A = £ r {0) is not 
zero, Cd+ T (ri) is the non-zero polynomial A(?i + 1) • • • (n + 
r), and c<j+ r (n — r) = \(n — (r — 1)) ■ ■ -n. Let then R = 
[0, . . . , r — 1] be the set of roots of the latter polynomial. 
For r < n < pd — 1, if (n mod p) is not in R, then u„ is the 
highest-index unknown appearing with a non-zero coefficient 
in C n -r\ we can then eliminate it, by expressing it in terms 
of the previous u m 's. 

The unknowns we cannot eliminate this way are u n , with 
n in 

A = [n \ < n < pd - 1, (nmodp) 6 R ]; 

the residual equations are C„_ r = 0, for n in B = Bi U B2, 
with 

Bi = [n\r<n<pd— 1 and (n mod p) £ R] 



1 < j < h we have 



and 



[ n I pd < n < (p + l)d + r - 1 ] . 



To determine the dimension of Q, and later on find a basis 
of it, we rewrite the residual equations using the residual 
unknowns. 

For n — ip + j in B\ , the unknowns present in C„- r are 
u ip+j ^ r ^ d , . . . ,u ip+ j. Of those, only u ip+ j_ r _ d , . . . , Uip_i 
need to be rewritten in terms of [u ?l | n € A]; the others 
already belong to this set. Thus, it is enough to express all 
iti p _, — d, . . . ,Ui p -i in terms of [u n | n G A], for 1 < i < d. 

For n in B2, the unknowns in C n -r are u n _, — d, ■ ■ • , Mpd-i 
(the higher index ones are zero). So, it is enough to compute 
u P d-r-d, ■ • • , u P d-i in terms of [u„ | n 6 A]. This is thus the 
same problem as above, for index i = d. 

Expressing all needed unknowns using A. Let A' = 

[0, . . . ,pd— 1] — A. For n in A, one can rewrite the equation 
C n -r = as the first order recurrence 



(5) 
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note that for n/0 mod p, A(n +p) = A(n). Let next B 
be the matrix factorial A(p — 1) • ■ • A(r). Then we have the 
equalities, for 1 < i < d: 
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U(i-l) p -d 








= B 
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Note that \A\ 


= dr; we 


let u 


3e the dr x 1 column-vector 



consisting of all u n , for n in A. Let further Co be the (d + 
r) x dr zero matrix. For 1 < i < d, suppose that we have 
determined (d + r) x dr matrices Ci, . . . , Ci_i such that, for 



C,u 



and 



with 



D, 



■J.jp-i 



Orxl, Ir 0,. 



Ujp-\-r — 1 



(6) 



ii = U ~ !) r and 4 = ( d "i " 1 ) r - 



Letting C' i _ 1 be the matrix made of the last d rows of Cj_i, 
we define 



Ci = B 



rxii. 



1 Ir 0rx<' 



then, (6) is satisfied at index i as well. 

Rewriting all residual equations using A. Combining 
all previous information, we obtain a matrix equality of the 
form u' = Du, where u' is the column vector with entries 

Ui p -, — d, , Uip+, — 1, for 1 < i < d, and where D is the 

matrix obtained by stacking up Di, . . . , D^. 

We have seen that all indeterminates appearing in the 
residual equations C„- r , with r in B, are actually in u'. By 
evaluating the coefficients Co, . . . , Cd+ r at n — r, for n in B, 
we obtain the matrix D' of the residual equations, expressed 
in terms of the unknowns in u'. Hence, the matrix E = D'D 
expresses the residual equations in terms of u n , for n in A. 

By construction, the dimension of Q equals the dimension 
of the nullspace of E. Knowing a basis of the nullspace of 
E, one deduces a basis of Q using (5), to compute all u n for 
n in A' . 

Cost analysis. By [6, Lemma 7], one can compute B in 



time Ti = 0(d"M(r i/2 p 1/ ^) log(rp)), which is 0(d t ^r 1/ V 7 )• 
Computing a matrix Ci requires one matrix multiplication 
of size (d + r,d) x (d, dr). In view of the inequality r < d, 
using block matrix multiplication, this can be done in time 
0(d" J r). Thus, computing all needed matrices Ci takes time 

T 2 = 0(d" +1 r)- 

The matrix D has size (d(d + 2r) x dr; no more com- 
putations are needed to fill its entries. The matrix D' has 
size d(r + 1) x d(d + 2r). Its entries are obtained by eval- 
uating Co, . . . , c r+ d at all n in B. Since deg(ci) < r and 
\B\ = d(r + 1), this takes time 0(M(dr) log(dr)) per polyno- 
mial. Since r < d, the total time is T 3 = 0(d 2 M(r) log(r)) G 
6(d 2 r). 

The matrix E = D'D has size d(r + 1) x dr; using block 
matrix multiplication with blocks of size dr, it can be com- 
puted in time T4 = 0(d" +1 r" _1 ). A basis of its nullspace 
can be computed in time T$ = 0(d u 'r u '). 

Given a vector [u n n G A] in the nullspace of E, one can 
reconstruct [u n < n < pd] using (5). This first requires 
evaluating all coefficients of all equations C n -r, for n in A' — 
[0, . . . ,pd - 1] - A, which takes time T 6 = 0(dM(s) log(s)), 
with s — max(r,pd). 

Then, for a given [u n \ n £ A] in the nullspace of E, de- 
ducing [u ?l < n < pd] requires \A'\ < pd matrix- vector 
products in size d + r. The dimension of the nullspace is 
O(dr); we process all vectors in the nullspace basis simulta- 
neously, so that we are left to do pd matrix products in size 
(d + r) x (d + r) by (d + r) x dr. The cost of each product 
is 0{d u} r), so the total cost is TV = 0(d" +1 rp). 

Summing Ti, . . . , Tg proves the first part of the proposi- 
tion. Adding to this T$ and T7 gives the second claim. 



3.3 Proof of Theorem 1 

Let T and Q be as above. By Proposition 1, &m\f p {Q) = 
if and only if dim F ^p) (T) = 0. Hence, the first estimate of 
Proposition 2 proves our first claim. 

Suppose that dmiF p (S) 7^ 0, and let ui, . . . , Uk be an F p - 
basis of Q. Proposition 1 implies that ui, ...,Uk generates 
T over ¥ p (x p ). We deduce an F p (a; p )-basis SB of T in a 
naive way: starting from SB = [ui], we successively try to 
add U2, ■ ■ ■ to SB. Independence tests are performed at each 
step, using the following lemma. 

Lemma 4. Given u\, . . . ,ui in ¥ p \x] of degree less than 
pd, one can determine whether they are linearly independent 
over ¥ p (x p ) in time 0(£" +2 dp). 

Proof. It suffices to compute their Wronskian determinant. 
The determinant of a matrix of size £ can be computed using 
0{f +1 ) sums and products [2]; since here all products can 
be truncated in degree idp, the cost is 0{f +1 U(£dp)). □ 

At all times, there are at most r elements in SB, so we always 
have £ < r + 1. Since we also have k < dr, the overall time 
is 0(d 2 r" +3 p), as claimed. 

4. P-CURVATURE: FIRST ORDER 

For first order operators, there is a closed form formula for 
the p-curvature. Let L = d — u, with u in ¥ p (x); then, by [25, 
Lemma 1.4.2], the p-curvature of L is the lxl matrix with 
entry u^ -1 ' + u p , where the first term is the derivative of 
order p — 1 of u. In this case, we do not distinguish between 
the p-curvature and its unique entry. 

The case of first order operators stands out as the only 
one where a cost polynomial in log(p) can be reached; this 
is possible since in this case, we only compute 0(d) non-zero 
coefficients. As per our convention, in the following state- 
ment, we take L not necessarily monic, but with polynomial 
coefficients. 

Theorem 2. Given L — ad — b in W p [x]{d) of bidegree 
[d, 1) that satisfies H, one can compute its p-curvature in 
time 0(dM(d) log(p)) C 0(d 2 log(p)). 

Proof. Since the p-curvature belongs to ¥ p (x p ), it suffices 
to compute its pth root. Computing the p-curvature itself 
requires no extra arithmetic operation, since taking p-powers 
is free over ¥ p , as far as arithmetic operations are concerned. 
Hence, we claim that the rational function 

{(a) + U ) ={(a) ) + a 

can be computed in time 0(dM(d) log(p)). Of course, the 
only non-trivial point is to compute s = (u^ p ~ 1 ^) 1 ^ p , with 
u = b/a. 

Observe that a p u ip - 1) is a polynomial of degree less than 
dp, so as is a polynomial of degree less than d. Hence, it 
is enough to compute the power series expansion s mod x d . 
From this, we deduce the polynomial as by a power series 
multiplication in degree d, and finally s by division by a. 

Let us write the power series expansion u = EooUiX 1 . 
Then, the series s equals — Ei>oUi p a; 8 , so it is enough to 
compute the coefficients (iii p )i<d. 

We start by computing the first coefficients uo, . . . , Ud-i 
by power series division, in time 0(M(d)). From these initial 



conditions, the coefficients u p , . . . , u p +d-i can be deduced 
for 0(M(d) log(p)) operations using binary powering tech- 
niques, see [16] or [3, Sect. 3.3.3]. Iterating this process d 
times, we obtain the values Ui P , . . . , Uip+d-i, for i < d, in 
time 0(dM(d)log(p)). □ 

As an aside, note that by Lemma 2, a rational function u is a 
logarithmic derivative in F p (a;) if and only if u^ -1 - 1 +u p = 0. 
This point also forms the basis of Niederreiter's algorithm 
for polynomial factoring [22]. 

5. P-CURVATURE: SECOND ORDER 

For second order operators, it is possible to exploit a certain 
linear differential system satisfied by the entries of the p- 
curvature matrix: already in [15, 26], one finds a third order 
linear differential equation satisfied by an anti-diagonal en- 
try of the p-curvature, for the case of operators of the form 
d 2 + s, or more generally d 2 + rd + s, when r^ -1 - 1 + r p = 0. 

In this section, we let L have the form vd 2 + wd + u, 
with u, v,w in ¥ p [x] of degree at most d. We assume that 
d > 2 and p > 2, and that H holds (we do not repeat these 
assumptions in the theorems); we let A be the companion 
matrix of L and let A p be its p-curvature. 

We give partial results regarding the computation of A p : 
we give algorithms of cost (5(p 1//2 ) or <D(p) to test proper- 
ties of A p , or compute it in some cases, up maybe to some 
indeterminacy. Though these algorithms do not solve all 
questions, they are still substantially faster than the ones 
for the general case in the next section. 

The trace of the p-curvature. We start by an easy but 
useful consequence of the result of the previous section: the 
trace of A p can be computed fast. 

Theorem 3. One can compute the trace r of A p in time 
0(dM(d) log(p)). 

Proof. The p-curvature of a determinant connection is the 
trace of the p-curvature of the original connection [21, 28]. 
Concretely, this means that the trace of A p is equal to the 
p-curvature of vd + w. By Theorem 2, it can be computed 
in time 0(d M(d) log(p)). □ 

Testing nilpotence. As a consequence of the previous 
theorems, we obtain a decision procedure for nilpotence. 

Corollary 1. One can decide whether A p is nilpotent 
m time 6{d"p l/2 + gT +1 ). 

Proof. The p-curvature A p is nilpotent if and only if its 
trace and determinant are both zero. By Theorem 3, the 
condition on the trace can be checked in time logarithmic 
in p. By Lemma 2, the second condition det(A p ) = is 
equivalent to the fact that Lu = has a non-zero solution, 
which can be tested in the requested time by Theorem 1. □ 

The eigenring. To state our further results, we need an 
extra object: the eigenring £(L) of L. This is the set of 
matrices B in ^2{¥ p (x)) that satisfy the matrix differential 
equation 

B'=BA- AB (7) 

(our definition differs slightly from the usual one in the 
sign convention). By construction, the eigenring £(L) is a 
F p (a; p )-vector space of dimension at most 4, which contains 
the p-curvature A p . Then, we let 7 be its dimension over 
¥ p (x p ); we will prove later on that 7 is in {2, 4}. 



Let further T be the set of solutions of Ly = in ¥ p {x) 
and let (3 be its dimension over ¥ p {x p ). Then, our main 
results are the following. 

Theorem 4. One can compute in time 0(d iJ+1 p) : 

1. the dimensions 7 G{2, 4} of£(L) and (3 €{0,1,2} of T ; 

2. A p , if 7 = 4 or /3 = 2. 

3. A p , up to a multiplicative constant in ¥ p [x p ] of degree 
at most pd, 1/7 = 2 and the trace t = 0. 

4- a list of two candidates for A p , if 7 = 2 and = 1. 

The rest of this section is devoted to prove this theorem. 

The dimension of the eigenring. The following lemmas 
restrict the possible dimension 7 of £(L). 

Lemma 5. If A p has the form AI2, then 7 = 4. 

Proof. In this case, the commutator of A p in ^2(¥ p (x)) is 
./^j(F p (x)) itself, so it has dimension 4 over F p (x). Then, [12, 
Prop. 3.5] implies that £(L) has dimension 4 over ¥ p (x p ). □ 

Lemma 6. Either 7 = 2, or 7 = 4. In the second case, 
A p is equal to J I2 , where t is the trace of A p . 

Proof. Corollary 1 of [12] shows that if the minimal and 
characteristic polynomials of A p coincide, then £(L) equals 
¥p(x p )[A p ]. In this case, ¥ p (x p )[A p ] has dimension 2 over 
¥ p (x p ). Else, the minimal polynomial of A p must have de- 
gree 1, so A p is necessarily equal to ^ I 2 , and we are under 
the assumptions of the previous lemma. □ 

Computing 7 and /3. The equality (7) gives a system of 
four linear differential equations of order one for the entries 
61,1 j • ■ • , 62,2 of B. An easy computation shows that (7) is 
equivalent to the system 

V 3 h2A + Ab' 2 A + B&2.1 = 0, 



v 2 b ia + Rb'li + S&a.i + Tb 2 ,i = 0, 



(8) 
(9) 



v(bl,l — 6 2 , 2 ) + V&2.1 — W&2.1 = 0, (10) 

61,1 + 62,2 = 0, (11) 
where A, B, R, S, T belong to F p [s], and are given by 
A — v(—2us'v + 2wv' + 4uv — w 2 ), 

B = vw{v" — w')+v'w(w— 2v')-\-2u'v 2 — 2vuv' — w"v 2 +2v'w'v 
and 

R = v 2 /2, S = -vw/2, T = v'w/2-vw'/2 + uv. 

Since Equation (11) is equivalent to 61,1 +62,2 £ ¥ p (x p ), we 
readily deduce that the dimension 7 of £(L) equals 7' + 1, 
where 7' is the dimension of the solution-set of (8). 

Computing both 7 and f3 can be done using Theorem 1, 
with respectively r = 3 or r = 2, and in degree respectively 
at most Ad or d. This proves point 1 of Theorem 4. 

If 7 = 4, we are in the second case of Lemma 6. Since the 
trace can be computed in time 0(d 2 log(p)) by Theorem 3, 
point 2 of Theorem 4 is established in this case. If fj = 2, 
then A p is zero by Lemma 2, so point 2 of Theorem 4 is 
established as well. 

Eigenrings of dimension 2. The rest of this section is 
devoted to analyze what happens if £(L) has dimension 



7 = 2 over ¥ p (x p ), so that the dimension 7' of the solution- 
space of (8) is 1. In this case, the information provided by 
the eigenring is not sufficient to completely determine the 
p-curvature. However, it is still possible to recover some 
useful partial information. To fix notation, we write the 
p-curvature as 

/i,i /l,2 
/2.1 /2,2 



Ap — 



Lemma 7. Ifj — 2,F = v p f2,i is a nonzero polynomial 
solution of degree at most pd of Equation (8) . 

Proof. Since the p-curvature A p belongs to the eigenring, 
its entries /i,i, . . . , / 2 ,2 satisfy (8) to (11). Lemma 1 shows 
that F = v p /2,i is a polynomial solution of degree at most pd 
of Equation (8). Moreover, F cannot be 0, since otherwise 
Equations (9) to (11) would imply that A p has the form AI 2 
for some A in ¥ p (x p ). By Lemma 5, this would contradict 
the assumption 7 = 2. □ 

Lemma 8. Suppose that 7 = 2 and let u £ ¥ p [x] be the 
nontrivial polynomial solution of minimal degree of Equa- 
tion (8). There exists a nonzero polynomial c in ¥ p [x p ] of 
degree at most pd, such that the entries of A p are given by 

^ = \( T+ VvG U - U '))' 
A' 2 = --pT5i Ru " + Su ' + Tu ) 1 

f c 

f2 - 2 = \{ t ~^G u ' u '))- 

Proof. By Lemma 7, the polynomials F and u both satisfy 
Equation (8); thus, they differ by an element c in ¥ p (x p ). 
Moreover, the minimality of the degree of u implies that 
c = F/u actually belongs to F p [o; p ] and has degree at most 
deg(F) < pd. The rest of the assertion follows from the 
relations F — v p f2,i, r = /i,i + /2, 2 and the equalities (9) 
and (10). □ 
Concluding the proof of Theorem 4. To conclude, we 
consider two special cases. If r = 0, as in [26], the pre- 
vious lemma shows that A p is known up to a multiplica- 
tive constant in F p [a; p ], as soon as the polynomial u has 
been computed. In this case, Corollary 1 shows that one 
can compute a non-zero solution Mo of (8) in the required 
time. The minimal degree solution u by clearing out the 
factor in F p [s: p ] in uo using [18, Ex. 14.27], in negligible 
time 0(M(dp) log(dp)) C O(dp), and the substitution in the 
former formulas takes time 0(dp) as well. 

If = 1, L has a non-trivial polynomial solution, so by 
Lemma 2 the determinant of A p is zero; the additional equa- 
tion /i,i/ 2l 2 = fi. 2/2,1, in conjunction with the formulas in 
Proposition 7, uniquely determines the polynomial c 2 and 
thus leaves us with only two possible candidates for A p . 

6. P-CURVATURE: HIGHER ORDER 

In this final section, we study operators of higher order, 
and we prove that the p-curvature can be computed in time 
subquadratic in p. 

Theorem 5. Given L in¥ p [x](d) of bidegree (d,r), one 
can compute its p-curvature in time 



u i2 2u/3 



0(r"d p 



+ r 



Hence, the exponent in p is 1 + < 1.79 < 2; in the best 
possible case ui = 2, we would obtain an exponent 5/3 in p, 
unfortunately still not optimal. 

As a result of independent interest, we also give an algo- 
rithm for computing the image of a matrix of rational func- 
tions by an differential operator similar in spirit to Brent 
and Kung's algorithm for modular composition [9]; to our 
knowledge, no prior non-trivial algorithm existed for this 
task. 

6.1 Preliminaries 

Euler's operator. Besides operators in the usual vari- 
ables x, d, it will also be convenient to consider operators in 
¥ p (x){9) or ^ n (¥ p (x)){9) , where 9 is Euler's operator xd, 
which satisfies the commutation rule 8x — x9 + x. To avoid 
confusion, we may say that L has bidegree (d, r) in d or in 
9, if L is written respectively on the bases (x, d) or (x,9). 

Conversion. Given an operator L in ^ n (¥ p [x]){d) of bide- 
gree (d, r), L' — x r L can be rewritten as an operator in 9 
with polynomial coefficients. The operator L' has bidegree 
(d + r,r) in 9. By [5, Section 3.3], computing the coefficients 
of L' takes time 0(n 2 (d + r)M(r) log(r)). Since representing 
all coefficients of L' requires 0(n 2 (d + r)r) elements, this is 
quasi-linear, up to logarithmic factors. 

Multiplication. Next, we give an algorithm for the multi- 
plication of operators with rational coefficients of a special 
type, inspired by that of [5] (which handles polynomial co- 
efficients). The algorithm relies on an evaluation / interpo- 
lation idea originally due to [24] , and introduces fast matrix 
multiplication to solve the problem. 

Lemma 9. Let b G ¥ p [x] be of degree at most d, with 
6(0) =fc and let j,p be in ¥ p (x)(d), with 

where gj andmj G ¥ p [x] have degrees at most d(h—j). Then 
if 2h < p — 1, one can compute n — 7/1 in time 0(h"d 2 ). 

Proof. Define n* = b 2h ?7,7* = b h -y and fj,* = b h /j,. A quick 
verification shows that these operators are in F p [x] (d) , of re- 
spective bidegrees bounded by (2dh, 2h), (dh, h) and (dh, h). 

We first compute j(x 3 ) and fi(x j ) mod x 2dh + h + 1 f or j < 
2h. This is done by computing the corresponding values of 
7* and p* , and dividing the results by b h . The former com- 
putation takes time 0(M(dh 2 )) using algorithm Eval of [5]; 
the latter 0(hM(dh)) by Newton iteration for power series 
division. Our assumption 2h < p — 1 ensures that divisions 
performed in the evaluation algorithm (and in the interpo- 
lation below) are well-defined. 

From the values of 7 and p, the values n(x 3 ) mod x 2dh+1 
are obtained as in [5, Th. 3]; the cost is 0(h"d 2 ). We can 
then compute the values of rf in time 0(hM(dh)) by fast 
polynomial multiplication. Knowing its values, we recover 
rj* using algorithm Interpol of [5]; this takes time 0(M(dh 2 )). 
Finally, we deduce n by division by b 2h ; this takes time 
0(hM(dh) log(dft)), using fast gcd computation. □ 

Left and right forms. Let L G ^ n (¥ p [x]){9) have the 
form 



with £i G ^ n (¥p[x]) of degrees at most d. It can be rewrit- 
ten 

L = t (x) + 9*£i{x) + ■■■ + 9 r t r (x), 

with £\ in ^ n (¥ p [x]) of degrees at most d as well. The 
former expression will be called the right-form of L; the 
latter is its left-form. 

Lemma 10. Let L have bidegree (d,r) in ^£ n {¥ p [x\)(9) , 
given in its right-form (resp. in its left-form). Then 
one can compute its left-form (resp. right-form) in time 
0(n 2 dM(r)log(r)) C 6(n 2 dr). 

Proof. We prove one direction only; the other is similar. 
Given the right-form of L, we can (without performing any 
operation) rewrite L = ^2j <d x 3 Li, where Lj has constant 
coefficients and order at most r. Since x 3 Li(9) — Li(9—j)x J , 
the result follows by using algorithms for polynomial shift 
by j [17]. □ 
The number of elements needed to represent L in either left- 
or right-form is 0(n 2 dr), so the previous algorithm is quasi- 
linear, up to logarithmic factors. 

6.2 Evaluation 

For L in JK r {¥ p [x]){d) or .^ r (¥ p [x}){9) and A in .# r (F p (x)), 
I/A denotes the matrix in ..4ff r (F p (x)) obtained by applying 
I to A. In this subsection, we give cost estimates on the 
computation of LA. 

The polynomial case. We start with the case of an opera- 
tor with polynomial coefficients, which we apply to a matrix 
with polynomial entries. We use an operator in 9, since this 
makes operations slightly more convenient than in d. As in 
Section 3, we make assumptions on the relative sizes of the 
input parameters (here S,p,e), for simplicity's sake. 

Lemma 11. Given L G ,£ r (¥ p [x])(9) of bidegree (S,p) 
and E G */£ r (¥ p [x\) of degree e, one can compute LE in time 
0(r" pe"- 2 8 A -^), assuming 5 G O(e) and e G 0(p 1/2 S). 

The cost can be rewritten as 0(r"pe (5/e) 3- "). Since lu < 3 
and 5 G 0(e), this is always better than O(r^pe): the cost 
ranges from 0(r"p5) for a hypothetical u> — 2 to 0(r <i, pe) 
for uj = 3. As a matter of comparison, let us write 

L = Y Jt < P ^e\ ^e.£r{¥ P [x\). 

Computing LE naively amounts to computing all (9 l E for 
i < p, multiplying them by the respective coefficients £i, and 
summing the results; the cost is in 0(r"pe), so our estimate 
is better. 

Proof. Our result is achieved using a baby steps / giant 
steps strategy inspired by Brent-Kung's algorithm for power 
series composition [9]. Let k = Lp 1 ^ 2 ] and h = \p/k~\. First, 
we rewrite L in left-form, as 

L = J2 i < p e i et(x); 

by Lemma 10, the cost is Ti =0(r 2 SM(p) log(p)) C 6{r 2 Sp). 
Next, L is cut into h slices of the form 

L +9 k L 1 + ---+9<- h - 1)k L h - 1 , i.e. I. y: ,0"l.,. 

Each Lj has order less than k and can be written as 



L = £o(x) + l±{x)9 + ■ ■ ■ + £ r (x)9 r , 



where for jk + i > p, £j k+i is zero. Finally, we rewrite each 
Lj in right-form: 

^■ = Ei< fc 4i( a, )* i ' ( 12 ) 

where all { have degree at most 5. By Lemma 10, the cost 
is T 2 = 0(hr 2 SM(k)\og(k)), which is in d(r 2 8p) as before. 
To apply L to E, we first compute the baby steps 

Eo = E, Ei = 6E, . . . , Efe_i = 9 k X E; 

then, we deduce all LjE, for j < h; finally, we do the giant 
steps 

All Ei can be computed in time T3 = 0(r 2 p 1//2 e), by suc- 
cessive applications of 9. The cost T4 of deducing the poly- 
nomials LjE is detailed below. Finally, one recovers LE 
by first computing all # jfc LjE, for j < h, and then sum- 
ming them. Since 1 (x-') — fx\ 6^ k can be applied to LjE 
in time 0(r 2 e log(p)), so the total cost of this final step is 
T 5 = 0{r 2 eh\og{p)) C 0{r 2 p 1/2 e). 

It remains to compute all LjE, given all Ei; we compute 
them all at once. In view of Equation (12), we have 

£j E = E I<fc ^j,i E i> 

where the E; are known. We cut Ei into slices of length 5: 

Ei = E„< s Ei, ua; *«, 

where Ei jU has degree less than 5 and s = \e/8] < 2e/5. 
This gives 

" E 4 E Ei,^ = E^ U E 

We will compute all inner sums 

5^i<fc ^i.jEi.u 

at once, for j < h and « < s; from this, one can recover all 
LjE in time 0(r 2 ep 1/2 ). 

The computation of these sums amounts to perform a 
(hx k) x (k x s) matrix multiplication, with entries that are 
polynomial matrices of size r and degree at most 8. Since 
e G 0(/9 1/2 <5), we have s G 0(p 1/2 ). Hence, we divide the 
previous matrices into blocks of size s and we are left to do 
a (0{p 1/2 /s) x 0(p 1/2 /s)) x (0(p 1/2 /s) x O(l)) product of 
such blocks, where p 1 ^ 2 /s is lower-bounded by a constant. 
Multiplying a single block takes time 0(r u s"M(6)), so the 
total time T 4 is 0(r u ps" ~ 2 M (<$)), which is 6(r>e a '- 2 <5 3 - a '). 

The conclusion of Lemma 11 comes after a few simplifi- 
cations, which shows that the dominant cost is T4, for the 
final linear algebra step. □ 

The rational function case. Next, we study the applica- 
tion of an operator to a matrix of rational functions A (we 
make some simplifying assumptions on the denominators 
in A, which will be satisfied in the cases in §6.3 where we ap- 
ply this result) . Besides, our operator is now in ^ r (¥ p [x] ) (d) 
rather than in ./£ r {W p [x\){8) . 

Because of the larger number of parameters appearing in 
the construction, the cost estimate unfortunately becomes 
more complex than in the polynomial case. 

Lemma 12. Let L G jV r {¥ v [x])(d) be of bidegree (S,p). 
Let A G .^ r (Fp(x)) be of the form B/b K , with b G W p [x] of 



degree at most d and B G ^ r (F p [x]) of degree at most nd. 
Define 

S' = S + p and e — (k + p)d + 5+1. 

7/6(0) 7^ and e G 0(p 1//2 <5'), one can compute LA in time 
6{r^pe^- 2 5' 3 ^). 

Proof. Let L' = x p L. Given L as an operator in d, we 
saw that we can write L' as an operator in 6, of bidegree 
(8' , p); the coefficients of L' in 6 can be computed in time 
0(r 2 SM(p) log(p)) C 0(r 2 5p). To conclude, it is enough to 
compute L'A, since then LA is deduced by a division by 
x p , which is free. 

For any i > 0, 6 Z A has the form Bi/fe K+l , with Bi in 
-# r (Fj, [x]) of degree at most (« + i)d. Thus, L'A has the 
form B*/b K+p , with B* G ^ r (¥ p [x]) of de gree less than e. 
with e = (k + p)d + 8' + 1. 

Knowing L'A mod x e , one can recover the numerator ma- 
trix B* through multiplication by b K+p ; a gcd computation 
finally gives L'A in normal form. These latter steps take 
time 0(r 2 M(e) log(e)) C 6(r 2 e). 

Since b(0) 7^ 0, the matrix E = A mod x £ is well-defined; 
it can be computed in time 0(r 2 M(e)) C (5(r 2 e) by power 
series division. Lemma 11 gives complexity estimates for 
computing L'E mod x E . Since this matrix coincides with 
L'A modulo x e , this concludes the proof of the lemma, 
as all previous costs are negligible compared to the one of 
Lemma 11. □ 

6.3 Computing the p-curvature 

Let L be in F p [x] {d) of bidegree (d, r) and let A be its com- 
panion matrix. We define the operator A G ^ r (F p (x))(d) 
as A = d + A; thus, as pointed out in the introduction, the 
p-curvature of L is obtained by applying A p_1 to A. 

To obtain a cost better than 0(p 2 ), we first compute a 
high enough power A' = A k of A; then, we apply A' to A 
k! times, with k! ~ (p — l)/k. Since p — 1 may not factor 
exactly as kk' , a few iterations of this process are needed. 

Computing A fc . Let £ = l r G F p [a;] be the leading co- 
efficient of L. Then, A has the form d + with A in 
.>£ r (F p [x]) and £ G F p [x] of degree at most d and £(0) 0. 
More generally, for k > 0, we can write A fc as 

A fc =E-=o A ^> with A ^ = ^J 
and £k.j in ^£ r {W p [x\) of degree at most d(k — j). 

Lemma 13. If k < p — 1, one can compute A k in time 
Oir^k^d 2 ). 

Proof. We use a divide-and-conquer scheme. Let h — 
[k/2\; we assume for simplicity that k = 2h; if k is odd, an 
extra (cheaper) multiplication by A is needed. We assume 
that A h is known, and we see it as an r x r matrix with en- 
tries that are scalar operators; hence, to compute A fe , we do 
O(r^) products of such scalar operators. All these products 
have the form 77 = 7p of the form seen in Lemma 9, so each 
of their costs is 0(h^d 2 ) = 0(k"d 2 ). □ 
. . / 

Computing A A. We fix k < p, and we compute the 
operators F — A k and T' = d k F. Writing k! = |_0 - l)/fcj, 
we compute the sequence 

A ( i)=A, A w =rA (< _ 1)) i = 2,...,k', 



so that A(fc/) = A kk A. Thus, we have A( fc /j = A k y, 
where the latter matrix is defined in Equation (1). Us- 
ing the subroutines seen before, a quick analysis not repro- 
duced here shows that the optimal choice is k = [(p— 1) 2//3 J . 
Then, computing F takes time 0(r"fc"d 2 ) = 0(r"p 2 " /3 d 2 ) 
by Lemma 13. 

By Lemma 1, each matrix A^) has the form B( i )/b lfe , 
with B(ij G ^# r (Fp[3;]) of degree at most dki. Given A^j, 
we compute A( i+i ) by first applying T to A(;) and dividing 
the result by d . 

The first step, applying V , is the more costly. We obtain 
its cost by applying Lemma 12, with (S, p) — (dk, k) and k = 
ik. Then, we have 5 G O(dk) and e G 0(idk). For all i < k , 
we are under the assumptions of that lemma; after a few 
simplifications, the cost becomes 0(r ul i u} ~ 2 dk 2 ). Summing 
over all i < k' , we obtain an overall cost of 0{r^k"^~ dk 2 ). 
Taking into account that k G 0(p 2/3 ) and k! G 0(p 1/3 ), this 
finally gives a cost of 0(r"dp 1+u ' /3 ) for computing A kk A. 

Computing the p-curvature. The definitions of k, k' im- 
ply that p - (p - 1) 2/3 < kk' < p - 1. To obtain the p- 
curvature A P_1 A, we iterate the previous process, replacing 
the required number of steps p — 1 by p — 1 — kk' , until the re- 
quired number of steps is O(l). Since p— 1 — kk' < (p— 1) 2//3 , 
it takes O(loglogp) iterations; hence, the overall time is still 
in 6{r"dp 1+ul/3 ). 
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